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Abstract 

The variance of the number of levels in an energy interval around a level with 
large quantum numbers (semiclassical quantization) is studied for a particle 
in a rectangular box. Sampling involves changing the ratio of the rectangle's 
sides while keeping the area constant. For sufficiently narrow intervals, one 
finds the usual linear growth with the width of the interval. For wider inter- 
vals, the variance undergoes large, non-decaying oscillations around what is 
expected to be the saturation value. These oscillations can be explained as a 
superposition of just a few harmonics that correspond to the shortest periodic 
orbits in the rectangle. The analytical and numerical results are in excellent 
agreement. 

I. INTRODUCTION 

Two decades ago, Casati, Chirikov and Guarneri [1] and Berry [2] have undertaken 
numerical and analytical studies of the level statistics in rectangular billiards for large quan- 
tum numbers (semiclassical quantization). Casati et al [1] demonstrated numerically that 
the distribution function of the nearest-neighbor energy level spacings is exponential (obeys 
Poisson law [3]), which is generally the case for classically integrable systems and is well 
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understood theoretically [4]. Further, they studied the behavior of the level rigidity [3] in 
the spectral staircase and showed that, following the initial linear growth with the width of 
the interval, the rigidity saturates to a constant value. This behavior was explained by Berry 
[2] who argued that the width of the interval at which the saturation occurs corresponds to 
the shortest periodic orbit in the billiard and that the saturation value of the rigidity can 
be obtained as a sum over the periodic orbits. Subsequently, the results of [1] have been 
reproduced with high numerical precision in Refs. [5] and [6]. 

In the previous work [7], we set out to interpret this result in terms of the global level 
rigidity and, in particular, proposed an ansatz for the level density correlation function 
that successfully describes the transition from linear to saturation behavior and that is also 
similar, aside from the energy scale, to its counterparts in the Gaussian ensembles (classically 
chaotic systems) [3], [4]. Furthermore, our correlation function led us to believe that the 
variance should exhibit interesting, albeit decaying, oscillations on approach to its saturation 
value (by the integral relationship [3] between the rigidity and the variance - see below - it 
should be twice that of the rigidity). 

Surprisingly, our numerical calculation indicated that the size of the oscillations of the 
variance does not decay with the increase of the width of the interval. Further investigation 
of the level correlation function revealed that a more accurate treatment of large energy 
scales gives an analytical explanation which is in excellent agreement with the numerical 
results. 

It should be emphasized that the key difference of our numerical procedure relative to [5] 
and [6] is that we averaged over the ratio of the rectangle's sides, as opposed to averaging over 
the spectrum keeping a fixed ratio of the sides. Indeed, both the onset of saturation and the 
saturation value scale as a square root of energy, that is depend critically on the position in 
the spectrum [2] . Consequently, averaging over the spectrum averages over these quantities. 
Notice, for instance, that the plot of variance in [5] does show oscillations. However, because 
of the averaging procedure, they appear to be decaying and of irregular pattern and did not 
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receive due acknowledgment. 



II. NOTATIONS AND KNOWN RESULTS 

We will consider intervals [e — E/2, e + E/2] , £ <C e, where the states with energies 
near e have large quantum numbers and can be described semiclassically. We will denote 
N (e) the spectral staircase [4] 

Af(e) = Y,0 (e-e*) (1) 

k 

where 9 is unit step function and k are energy eigenstates. More precisely, we will consider 
the 2D-reduced [4] spectrum where the systematic energy dependence of Af (e) has been 
eliminated 2 . Towards this end, we make the substitution 

e-CA/-( e )> (2) 

and, in particular, — > (Af (£&)), where (Af (s)) is the ensemble average of N '(e) (in the 
limit of an infinitely large ensemble). With such definition of energy (reduced spectrum), we 
have 

(Af(e))^(Af(e))=e = (p(e))e (3) 

that is (p(e)) = 1, where 

k 



x We learned about Ref. [5] after completion of this work. 

2 For a finite domain, in particular, we must eliminate the systematic energy dependence of (M (e)) 
due to boundary corrections [4], [8], [9]; for a rectangle this dependence can be obtained directly 
[9], [10] from the energy levels of a particle in a box, eq. (7) below, and is given by (AT (e)} = 
e - 2(3^JeJn + 1/4, where j3 = (a 1 / 4 + a~ 1 ^) /2, a is defined in eq. (8) below and e is measured 
in units of A = 2Trh 2 /mA, eq. (15) below. 
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is the density of states. In dimensional units, (3) can be written as 
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where 



(P(e)) 



mA 



(6) 



m/2iih 2 being the 2D density of states in the thermodynamic limit and A the system area. 
For a rectangle, the energy eigenvalues are 



and averaging is performed over the "a-ensemble," that is by sampling the values of 

a = pr ( § ) 

where a and b are the rectangle's sides, while keeping A = ab, and consequently (p(e)), 
constant. Numerically, we use algebraic numbers for a to avoid accidental level degeneracies 
(we don't use transcendental numbers because they are too close to rational [4]). 

In this work we will not discuss the distribution of the nearest-neighbor energy level 
spacings, except to state that our computations [10] are congruent with those of [1], [5] 
and [6]. Here, we will study the correlations and fluctuations in the energy spectrum and, 
towards this end, the following measures will be used. First, it is the "level rigidity" A 3 , 
defined [2], [3] as the best linear fit of the spectral staircase on the interval [e — E/2,e + E/2] 




(7) 





Fig. 1 shows an example of such a fit. 



For the number of levels iV on the interval [e — E/2, e + E/2] 




(10) 



the variance 
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E(e,£) = <(iV-(iV» 2 ) (11) 
is another measure of the fluctuations. Clearly, 

(AT) = (p) E (12) 

for an infinitely large ensemble. 

Introducing the mean level spacing (A) as 

m = ^{ji) (13) 

we find that, due to the iV-fluctuations S oc E over the intervals E narrower than a certain 
energy scale, the following relationship holds [10] 



(N) + °{(NY 



<P>(A> = 1 + — + 0(— 1 ) (14) 



For classically integrable billiards the latter holds for E <C \J e (p) 1 (see below), so that 
we can neglect the difference between (A) and only for sufficiently large intervals 

(N) 3> 1, in which case we can introduce the limiting value of mean level spacing ("global 
mean") as 

A = JL = H^ (15, 
(p) mA 

In what follows we will, unless mentioned otherwise, measure energies in units of A (wherein 
A = (p) 1 = 1, as in eq. (3)). 

General relationships can be established between the fluctuations measures and the cor- 
relation function of the level density [3], 

K(e 1 ,e 2 ) =(5p(e 1 )5p(e 2 )) (16) 
8p(e) =p(e)-(p(s)) (17) 

regardless of the specific form of K (e±, £2), for instance, 

re+E/2 re+E/2 

E{e,E)= / K(e 1 ,e 2 )de 1 de 2 (18) 

Je-E/2 Je-E/2 



Using these relationships one can further show that £ supersedes A 3 via an integral rela- 
tionship [3] 

A 3 (e, = J ' dx (E 3 - 2xE 2 + x 2 ) £ (e, x) (19) 

which, again, does not depend on the particulars of the level statistics and level correlations 
in the semiclassical energy spectrum. 

For Gaussian ensembles, corresponding to classically chaotic ergodic systems, the func- 
tional form of the level correlation function K (si, e 2 ) is well understood [3]. Denoting 

e = — - — , v = e 2 - e\ (20) 

we find that it can be written, in most general terms, as 

K (e 1 ,e 2 ) =K(u) = A~ 2 [d (|) - K (|)] (21) 

where 

/•oo 
K (x) dx = 1 (22) 
-oo 

and 

/oo 
K (u) dcu = (23) 
-oo 

The interpretation of these formulas are as follows. The 5-function term in (21) describes 
uncorrelated levels. The /C-function describes "level repulsion" and has a scale of A. The 
integrals (22) and (23) converge to their values on the scale of uo ~ A and reflect the fact 
that an overall "level rigidity" develops on the scale of A. 
Indeed, by definition of 5p (e), 

J 5p (e) de = (24) 

and thus 

J K ( £l , e 2 ) dei = J K ( £l , e 2 ) de 2 = (25) 



where integration is over the entire energy spectrum. The integral can be converted to the 
form (23) when, given a position in the spectrum e, the scale on which the level rigidity (24) 
develops is much smaller than e. Eqs. (21)-(23) show that Gaussian ensembles become rigid 
on the scale of A. This fact can be also observed from the behavior of E 

£(£,£) =|,£<A (26) 
X(e,E) = Cln(£),£»A (27) 

where the constant C depends on the specifics of a Gaussian ensemble [3]. The behavior of 
A 3 , 

A 3 (£,£) = ^ E ^ A ( 28 ) 
A 3 (e,E) =^ In (£),£» A (29) 

is then easily recovered from that of £ using (19). Notice that the linear terms in eqs. (26) 
and (28) originate in the 5-function term in (21) (uncorrelated levels). 

In the semiclassical approach, Sp(e) is an oscillatory term that can be written as a sum 
over all periodic orbits [4]. Berry [2] used this approach to obtain the following limiting 
behaviors of A 3 (e, E) for the rectangular box: 

A 3 (e,E) = 1|, E < E max = (na^eAf 2 (30) 

A 3 (e,£) =4t?(x) 1/2 E E 3«^>^U (31) 

-> 0.0947y^, a ~ 1 (32) 

where the second term gives the "saturation rigidity" A^ and 

if Mi = M 2 = 

5m = 1/4 if one of M x and M 2 is zero (33) 

1 otherwise 

Here Mi and M 2 are the "winding numbers" of the classical periodic orbits [2] and the 
factor 5m differentiates between the self-retracing and non-self-retracing orbits. Notice that 



both Af and E max are the functions of the position e in the spectrum (oc s/s), which is in 
stark contrast with Gaussian ensembles. The quantum scale E max for the onset of saturation 
corresponds to the time of traversal of the shortest classical periodic orbit [2], 

£max = (34) 
-L min 

whose length is just twice the length of the rectangle's smaller side and the winding numbers 
are M x = 0, M 2 = 1. 

Notice that saturation of A 3 to Af in a rectangle, eqs.(30) and (31), is analogous to 
"saturation" of A 3 to the weak logarithmic dependence for a Gaussian ensemble, eqs. (28) 
and (29). Consequently, one should expect the level density correlation function for a rect- 
angular box to have the form similar to (21), except that the scale of /C is set by ~ y/eA. 
Indeed, in our previous work [7] we introduced a simple ansatz that resulted in the following 
expression for K: 



K(e,u) = A" 



\A/ 7TUJ \E max J 



(35) 



In this form, K (e, to) satisfies eq. (23) and is, thus, consistent with the level rigidity develop- 
ing on the scale of E max . Evaluation of A 3 using eq. (35) successfully reproduces saturation 
to A3 on this scale also. For E, eq. (35) predicts [7] eventual saturation to S°° = 2Af , 
but in an oscillatory fashion 

= sin (E/E max ) 
E/E mSuX 

(The oscillatory term does not contribute to A^ in eq. (19)). 

While successfully describing the onset of the level rigidity and the oscillations of £ on 
the scales ~ E max , the ansatz (35) does not capture the correct behavior of E on scales 
> E mSLX . The central result of this work is that the magnitude of the oscillations around 
S°° is not decaying and that the onset of the level rigidity on the scale of E m3iX is not, 
in fact, precise. This will be shown both numerically and analytically, the latter being 
based on a more careful analysis of the level correlation function obtained, alternatively, 
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either semiclassically or by the direct use of the explicit form of the level spectrum in the 
rectangular box. 



III. CORRELATION FUNCTION OF THE LEVEL DENSITY 

The analytical expression for the level density correlation function, in the form of the 
infinite sum, can be obtained in two ways. First, it can be recovered from the semiclassical 
derivation of Berry in Ref. [2], as per eqs. (23), (38), (58), and (60) there. Indeed, taking 
the inverse Fourier transform in the last of these equations, we find 



oo oo 



K{e,u) = ; y y 45m VV '- (37) 



1 



oo oo 



4 _ _ cos (yg {Ml + Mf y) 2 ~ cos (y^Mc) 



Mi>0M 2 >0 V ^"1 1 ±v± 2 M>0 



J Q~l 

The paragraph following eq. (38) in [2] explains the nature of the factor A5 M in the above 



(38) 



equation; our interpretation is that the factor of 4 is the intensity factor due to constructive 
interference between the closed orbit and its time-reversed for non-self-retracing orbits (both 
Mi and M 2 are non-zero). The latter is crucial for understanding of the change in level 
correlations due to time- reversal symmetry breaking (for instance, by a magnetic field). 

Alternatively, one can obtain (37) starting directly from eqs. (4) and (7). In Ref. [11], 
such a formalism was developed for 5p and was generalized to K (e, cu) in Ref. [7]. It is based 
on the use of the Poisson summation formula, whose net effect is to convert the sum over 
actual levels to the sum over classical periodic orbits. Formally, it is accomplished [11], [7] 
by converting the sum on (n, m) to the sum on (Mi, M 2 ) of the Fourier transforms and by 
neglecting the rapidly oscillating terms, that is harmonic terms whose arguments depend 
on e and change by ~ a/s/A 3> 1 when Mi j2 change by 1. As a result one recovers 3 [7] eq. 



3 Notice that the correct power in the denominator of eq. (20) in Ref. [7] is 1, not 1/2. Also, we 
are presently investigating if the coefficient 5m is properly recovered with this technique. 
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(37). 

Notice that in the limit of uo <C VeA, the summation over (Mi, M 2 ) can be replaced by 
integration over (x, y) in the entire plane, where 

x = M^^/^A y = M2a - 1/4 yM (39 ) 
Converting to polar coordinates, we find 

k ^ ^ ^ r "* sy^'h^iid < 4 °> 

that is the first term in eq. (35). Notice, however, that in reality the lower limit of p- 
integration should not extend to zero since Mi )2 are not simultaneously zero. The second 
(repulsion) term in eq. (35) is then obtained using the ansatz [7] in which the lower cut-off 
in f°° is chosen at p min = y min = a" 1/4 a/IttA/^. 

In the opposite limit u <C VsA, the terms in the series (37) become rapidly oscillating 
with the amplitude rapidly decreasing with the increase of the winding numbers. Conse- 
quently, the sum will be dominated by a just few terms with the smallest M 12 . The latter 
is the reason behind the non-decaying oscillations of E that will be discussed below. 



IV. LEVEL NUMBER VARIANCE 



Combining eqs. (18) and (37), we find 



e 

TT^A 



~ » sin 2 (EyJ% (M* a i/* + M*crW)) 

2^ 2^ 4S m— 



Mi=0 M 2 =0 



(M 2 aV2 + M|a-V2) 3 / 2 



7T 



5 A 



- - sin 2 [E^ (M 2 + Mi 



sin 2 {E^/^M) 



M>0 



M 3 



(41) 



(42) 



J a~l 



For narrow intervals, eq. (41) reduces to 



E 



(43) 



which follows either from eqs. (18) and (40) or directly from (41) if summation is replaced 
with integration using variables (39). It is also consistent with eqs. (19) and (30). 



10 



Fig. 2 shows the result of evaluation of E using eq. (41) (with the upper limit of 
summation on M 12 limited to 100) versus the numerical evaluation of E for the ensemble of 
200 algebraic a G [1,2]. We also show E obtained using the ansatz (35) of Ref. [7]. Clearly, 
the latter adequately describes the transition from the small-scale E <C E max linear behavior 
(43) to scales of E ~ E mSLX , but fails to describe the large-scale behavior, E 3> E max . On the 
other hand, the agreement of numerical evaluation with the theoretical result (41) for all E 
is quite remarkable. 

Two comments are in order. First, it was necessary to have a sufficient spread of a- 
values (between 1 and 2 here) to obtain reliable statistics for large quantum numbers with 
relatively small a-ensembles (200 values here). However, the dependence on a of the relevant 
parameters is quite weak, through c^ 1 / 4 , so for a E [1, 2] the difference from an ensemble 
whose a- values would be close to a fixed value (1 in this particular case) is not significant. 
Second, for E 3> -Emax it is sufficient to limit summation on M ij2 to 2 and 3 in the double and 
single sum respectively to obtain a curve which is very close to the full theoretical curve; as 
was already explained above this is because the amplitudes of the quickly oscillating terms 
in the sum rapidly fall off with the increase of M lj2 . 

We now turn to the theoretical curve in Fig. 3 that represents the level number variance 
given by (42). It is intended to describe the a-ensemble such that a ~ 1 and it can be 
effectively reproduced by the superposition of only 6 harmonics: 3 from the double some and 
3 (with commensurate frequencies) from the single sum. These are non-decaying oscillations 
whose amplitude is given by 

y± o 

=r~ ~ = (44) 
E°° 3 v ' 

where E + = max{E} — E°°, E~ = E°° — min{E} and E°° is the mean value of oscillating 
E which is obtained by substituting sin 2 — > 1/2 in eq. (41). The small asymmetry of the 
positive and negative swings of E-oscillations, E + ^ E _ , is due to the asymmetry of the 
Clausen function [12] Cls(x), 
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Cl 3 (x) = J2 



oo , 

cos kx 



k 3 



k=i 

for which max{C/3(x)} / min {Cl^(x)} m 1.3. This function originates in the single sum in 
eq. (42), which corresponds to self-retracing orbits. 

It is important to point out that if one averages over the range of a-values, or over the 
range of e- values [5], the beats that result from a superposition of a continuous range of 
harmonics will present themselves, over relevant interval widths, as decaying oscillations. 
This is due to the corresponding continuous range of the beat frequencies represented by 
their envelopes. 

We now turn to breaking of time reversal symmetry for charged particles due to a mag- 
netic field. Clearly, the condition for the latter is given by 

BA ~ 0o (45) 

where O — hc/e is the flux quantum. For such fields, the Larmor radius is much greater 
than either side of the rectangle 

„ mcv /— [~e . . 

r =7b~ vi U (46) 

where mv 2 /2 = e. Therefore, the deviation from the free specular scattering will be small 
and we can find both the level correlation function and the level number variance via simple 
substitution 

if Mi = M 2 = 
5m 5' M = 1/4 if one of M 1 and M 2 is zero (47) 
1/2 otherwise 

in eqs. (37) and (41). For the a ~ 1-ensemble considered above, E' is shown in Fig. 4. In 
this case 



Voo 



E 



oc 



0.7 (48) 

e' f - _ e;- 

E 7 



- 0.77, — — w 0.64 (49) 
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The increased asymmetry of E-oscillations underscores the increased relative contribution 
of the single sum to E in eq. (42) (self-retracing orbits). 

V. DISCUSSION 

The central results of this work are summarized by Figs. 2-4. The first of these graphs 
shows that the numerical evaluation of the level number variance for a particle in a rect- 
angular box is in excellent agreement with the theoretical result given by eq. (41). The 
second indicates a non-decaying oscillatory behavior of E. The third shows E when the 
time-reversal symmetry is broken. Figs. 3 and 4 can be successfully reproduced with just a 
small number of lowest harmonics in (42). 

It is remarkable that, given the center of the interval e, E (e, E) exhibits large, repro- 
ducible oscillations as a function of the interval width E. The main implication of this result 
is that while the level rigidity indeed develops on the scale set by ~ v^A <C s, it is accurate 
only in an approximation where a harmonic is replaced by its average - zero. In other words, 
(23) converges to zero on the scale of u ~ v^A only up to a sum of harmonic terms, as 
implied by eq. (37). 

We point out that E-oscillations are entirely consistent with the near straight line sat- 
uration of A 3 . In fact, A 3 also exhibit oscillatory behavior around A3 , but it is both 
parametrically small in A/e and decays as a function of E. This reduction of A3 can be 
traced, for instance, to integration in eq. (19). 

In the future, we will investigate the level number variance in a variety of finite-domain 
and potential problems. We will also address the explicit dependence of the level correla- 
tion function on the magnetic field for a detailed description of the time-reversal breaking 
transition and to address the orbital magnetism of non-resonant integrable systems. 
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VII. FIGURE CAPTIONS 



1. Figure 1 

Level rigidity A 3 (e,E) vs. the interval width E for e/A — 10 4 and 200 algebraic a G 
[1,2]. The straight line is (E/15A) and the two horizontal lines correspond to Berry's 
saturation rigidity for a = 1 (lower line) and a = 2. Notice that our saturation rigidity 
is slightly below the former, while expected to be between the two, and that the difference 
from is amplified by the factor of -\fsfK. 

2. Figure 2 

Numerical evaluation of £ (e, E) (red line), theoretical evaluation per eq. (41) (green 
line), and theoretical evaluation based on ansatz (35) (blue line) for e/A = 10 4 and 200 
algebraic a G [1, 2]. 

3. Figure 3 

£ (e, E) per eq. (42) for e/A = 10 5 . 

4- Figure 4 

T,(e,E) per eqs. (42) and (47) for e/A = 10 5 for the case of broken time reversal 
symmetry. 
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